####################################################################
# Replication code for
#   Lyall, Shiraito, and Imai, "Coethnic Bias and Wartime Informing"
#
# Code to create Table 7 in the Supplemental Appendix
#
# Author: Yuki Shiraito
# Created: March 2, 2015
####################################################################

rm(list = ls())

library(coda)
library(xtable)

load("misc_retaliation.RData")
load("Posterior_oprobit.RData")

mcmc.colnames <- colnames(results.saved[[1]][[1]])
sample.var.params <- as.mcmc(rbind(results.saved[[1]][[1]][, c("sigma2", "gamma2", "gamma3", "gamma4")],
                                   results.saved[[1]][[2]][, c("sigma2", "gamma2", "gamma3", "gamma4")],
                                   results.saved[[1]][[3]][, c("sigma2", "gamma2", "gamma3", "gamma4")]))
sample.indiv.coef <- as.mcmc(rbind(results.saved[[1]][[1]][, mcmc.colnames %in% varname.indiv],
                                   results.saved[[1]][[2]][, mcmc.colnames %in% varname.indiv],
                                   results.saved[[1]][[3]][, mcmc.colnames %in% varname.indiv]))
sample.vil.coef <- as.mcmc(rbind(results.saved[[1]][[1]][, mcmc.colnames %in% varname.village],
                                 results.saved[[1]][[2]][, mcmc.colnames %in% varname.village],
                                 results.saved[[1]][[3]][, mcmc.colnames %in% varname.village]))

summary.village <- summary(sample.vil.coef)
summary.var.params <- summary(sample.var.params)

col.pashtun <- grep("endorserPashtun", mcmc.colnames)
col.tajik <- grep("endorserTajik", mcmc.colnames)
indiv.coef.ff <- sample.indiv.coef[, -c(col.pashtun, col.tajik)]
summary.ff <- summary(indiv.coef.ff)
indiv.coef.pashtun <- sample.indiv.coef[, col.pashtun]
indiv.coef.pashtun[, 2:ncol(indiv.coef.pashtun)] <- indiv.coef.pashtun[, 2:ncol(indiv.coef.pashtun)] +
    indiv.coef.ff
summary.pashtun <- summary(indiv.coef.pashtun)
indiv.coef.tajik <- sample.indiv.coef[, col.tajik]
indiv.coef.tajik[, 2:ncol(indiv.coef.tajik)] <- indiv.coef.pashtun[, 2:ncol(indiv.coef.tajik)] +
    indiv.coef.ff
summary.tajik <- summary(indiv.coef.tajik)

table.vil <- cbind(summary.village$statistics[, "Mean"],
                   summary.village$quantiles[, c("2.5%", "97.5%")])
table.vars <- cbind(summary.var.params$statistics[,"Mean"],
                    summary.var.params$quantiles[, c("2.5%", "97.5%")])
table.ff <- cbind(summary.ff$statistics[, "Mean"],
                  summary.ff$quantiles[, c("2.5%", "97.5%")])
table.pashtun <- cbind(summary.pashtun$statistics[, "Mean"],
                       summary.pashtun$quantiles[, c("2.5%", "97.5%")])
table.tajik <- cbind(summary.tajik$statistics[, "Mean"],
                     summary.tajik$quantiles[, c("2.5%", "97.5%")])

rownames(table.vars) <- c("Variance of Village Random Effects",
                          "Cutpoint (``somewhat unlikely\'\' - ``neither unlikely nor likely\'\')",
                          "Cutpoint (``neither unlikely nor likely\'\' - ``somewhat likely\'\')",
                          "Cutpoint (``somewhat likely\'\' - ``highly likely\'\')")
colnames(table.vars) <- c("Mean", "2.5%", "97.5%")
xtable(table.vars)

baseline.table <- rbind(table.ff, table.vil)
colnames(baseline.table) <- c("Mean", "2.5%", "97.5%")
rownames(baseline.table) <- c("Other Pashtun", "Tajik", "Monthly income",
                              "Age", "Years in school", "Daily contact with ISAF",
                              "Daily contact with ANP",
                              "Daily contact with Other Ethnics",

                              "Property damage by ISAF (Pro-Taliban Pashtun)",
                              "Harm by ISAF (Pro-Taliban Pashtun)",
                              "Both damage by ISAF (Pro-Taliban Pashtun)",

                              "Property damage by ISAF (Other Pashtun)",
                              "Harm by ISAF (Other Pashtun)",
                              "Both damage by ISAF (Other Pashtun)",

                              "Property damage by ISAF (Tajik)",
                              "Harm by ISAF (Tajik)",
                              "Both damage by ISAF (Tajik)",
                              
                              "Property damage by ANP (Pro-Taliban Pashtun)",
                              "Harm by ANP (Pro-Taliban Pashtun) ",
                              "Both damage by ANP (Pro-Taliban Pashtun)",

                              "Property damage by ANP (Other Pashtun)",
                              "Harm by ANP (Other Pashtun) ",
                              "Both damage by ANP (Other Pashtun)",

                              "Property damage by ANP (Tajik)",
                              "Harm by ANP (Tajik) ",
                              "Both damage by ANP (Tajik)",

                              "Harm by Taliban (Pro-Taliban Pashtun)",
                              "Both damage by Taliban (Pro-Taliban Pashtun)",

                              "Harm by Taliban (Other Pashtun)",
                              "Both damage by Taliban (Other Pashtun)",

                              "Property damage by Taliban (Tajik)",
                              "Harm by Taliban (Tajik)",
                              "Both damage by Taliban (Tajik)",

                              "Other Pashtun times Contact (ISAF)",
                              "Tajik times Contact (ISAF)",

                              "Other Pashtun times Contact (ANP)",
                              "Tajik times Contact (ANP)",
                              "Intercept",
                              "SIGACTS by ISAF within 2km",
                              "SIGACTS by Taliban within 2km",
                              "iMMAP by Unknown within 2km",
                              "iMMAP by ISAF within 2km",
                              "iMMAP by Taliban within 2km",
                              "NSP aid sites within village",
                              "Military Installations within 2km",
                              "Altitude (log)", "Population (log)",
                              "% of Pashtun", "Control by Taliban", "Contested")
xtable(baseline.table)

endorsers.table <- cbind(table.pashtun, table.tajik)
rownames(endorsers.table) <- c("Endorser (Pashtun or Tajik)",
                               "Other Pashtun", "Tajik", "Monthly income",
                               "Age", "Years in school", "Daily contact with ISAF",
                               "Daily contact with ANP",
                               "Daily contact with Other Ethnics",

                               "Property damage by ISAF (Pro-Taliban Pashtun)",
                               "Harm by ISAF (Pro-Taliban Pashtun)",
                               "Both damage by ISAF (Pro-Taliban Pashtun)",

                               "Property damage by ISAF (Other Pashtun)",
                               "Harm by ISAF (Other Pashtun)",
                               "Both damage by ISAF (Other Pashtun)",

                               "Property damage by ISAF (Tajik)",
                               "Harm by ISAF (Tajik)",
                               "Both damage by ISAF (Tajik)",
                               
                               "Property damage by ANP (Pro-Taliban Pashtun)",
                               "Harm by ANP (Pro-Taliban Pashtun) ",
                               "Both damage by ANP (Pro-Taliban Pashtun)",

                               "Property damage by ANP (Other Pashtun)",
                               "Harm by ANP (Other Pashtun) ",
                               "Both damage by ANP (Other Pashtun)",

                               "Property damage by ANP (Tajik)",
                               "Harm by ANP (Tajik) ",
                               "Both damage by ANP (Tajik)",

                               "Harm by Taliban (Pro-Taliban Pashtun)",
                               "Both damage by Taliban (Pro-Taliban Pashtun)",

                               "Harm by Taliban (Other Pashtun)",
                               "Both damage by Taliban (Other Pashtun)",

                               "Property damage by Taliban (Tajik)",
                               "Harm by Taliban (Tajik)",
                               "Both damage by Taliban (Tajik)",

                               "Other Pashtun times Contact (ISAF)",
                               "Tajik times Contact (ISAF)",

                               "Other Pashtun times Contact (ANP)",
                               "Tajik times Contact (ANP)")
colnames(endorsers.table) <- rep(c("Mean", "2.5%", "97.5%"), 2)
xtable(endorsers.table)
